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The proliferation and migration dichotomy of the tumor cell invasion is examined within a two- 
component continuous time random walk (CTRW) model. The balance equations for the cancer 
cells of two phenotypes with random switching between cell proliferation and migration are derived. 
The transport of tumor cells is formulated in terms of the CTRW with an arbitrary waiting time 
QQ , distribution law, while proliferation is modelled by a logistic growth. The overall rate of tumor cell 

• invasion for normal diffusion and subdiffusion is determined. 
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\ One of the main features of malignant brain cancer is the ability of tumor cells to invade the normal tissue away 
from the multi-cell tumor core. Invasion of healthy tissue by a solid tumor (the core), and the role of oxygen and 
I— 1[ nutrient delivery have been the subject of extensive studies reflected in modern surveys, see e.g., [1]. Experimental 
t data for a glioma cancer show that the proliferation rate of migratory cells is lower in the invasion region than in the 
' core. It turns out that the proliferation and migration of cells are mutually exclusive: the high motility suppresses cell 
• . proliferation and vice versa. This phenomenon is known as the migration-proliferation dichotomy [2, 3]. The exact 
. ' mechanism of switching between the two phenotypes of glioma cells is not known. There are several phenomenological 
models for this dichotomy. One can assume that the diffusion coefficient of cancer cells is a decreasing function of 



> 
o 



PACS numbers: 05.40.-a, 05.40.Fb, 87.15.Vv, 87.17.Ee, 82.39.Rt 



I. INTRODUCTION 



cell density [4] . As a result the cancer cell motility is greater in the invasion zone because of the small density of cells 
there. One can also assume the dependence of the proliferation term on cell density such that the proliferation rate 
increases with density [5]. An interesting dynamical model for the phenotype switch was suggested in [6]. However, 
this mathematical model involves many parameters, some of which are difficult to estimate. Recently the authors 
proposed a stochastic approach for the proliferation-migration switching that involves only two parameters [7]. The 
transport process was formulated in the framework of the continuous time random walk (CTRW) [8-10]. The main 
, reason for employing the CTRW model was to give the mesoscopic description of cancer cell motility in terms of 
the random jump distribution and waiting times. One of the main purposes was to take into account anomalous 
transport (subdiffusion) leading to slow motility of cancer cells in the invasive zone. Among all possible cancer cell 
genotypes, leading to six main alternations of malignant growth [11], cell motility and invasion are most important for 
J — I our consideration. The standard diffusion approximation for the transport (which is the parabolic limit of kinetics) 
, together with a logistic growth yields an overestimation of the overall growth [12, 13]. Since the motility is the most 
critical feature of brain cancer, causing treatment failure, there is a need for a proper description of cancer cell motility 
beyond the standard diffusion approximation. In this connection, the hyperbolic limit of the multi-cellular microscopic 
' system is important [14] to take into account cellular interaction in the description of macroscopic dynamics. A very 
JH [ interesting agent-based model was developed recently by Mansury and Deisboeck [15]. The transport process is 
described in terms of the local-search mechanism performed by tumor cells. The purpose of this "conscious" search 
is to find and then invade the most permissive location in extracellular matrix. A simplified scheme of migration- 
proliferation dichotomy in terms of CTRW was considered in [10, 16]. It involves two steps: cell fission with the 
characteristic time 7} and cell transport with duration 7^. During the time scale 7}, the cells interact strongly and 
motility of the cells is small. During the time %, interaction between the cells is weak and motility of the cells is 
determined by a "jump" length ~ %. 

Cell invasion is a very complex process controlled by matrix adhesion (see review [2]). It involves several steps 
including receptor-mediated adhesion of cells to extracellular matrix (ECM), matrix degradation by tumor-secreted 
proteases (proteolysis), detachment from ECM adhesion sites, and active invasion into intercellular space created by 
protease degradation. One of the purposes of this paper is to give a description of this complicated cell transport 
in terms of a non-symmetrical random walk model with memory effects. Chemotaxis and haptotaxis are taken into 
account by the biased random walk of cells that respond to external signals without alteration and migrate away from 
the tumor core. Matrix adhesion effects are modelled by using the heavy-tailed waiting time distributions that lead 
to subdiffusion of tumor cells. 
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II. TWO-COMPONENT CTRW WITH PROLIFERATION 
A. Balance equations 

In this paper we present a detailed analysis of the migration and proliferation of glioma cells in the framework of a 
two-component continuous time random walk with proliferation. The paper is an essential extension of our Letter [7] 

with new results and examples. Based on experimental observations of migration-proliferation dichotomy, wc assume 
that the process of tumor cell invasion consists of two states. In state 1 (migratory phenotype) the cells randomly 
move but there is no cell proliferation. In state 2 (proliferating phenotype) the cancer cells do not migrate and only 
proliferation takes place. To describe the random switching between the two phenotypes, we employ the two-state 
Markov chain model. The cell of type 1 remains in state 1 during a waiting time ri and then switches to a cell of 
type 2. After a waiting time T2, spent in state 2, it switches back to a cell of type 1. Both waiting times n and T2 
are mutually independent random variables exponentially distributed with parameters /3i and /32: 

V{Tk) = l3kexp{-l3kTk) k = l,2. (1) 

Here the parameters (3k are the switching rates, namely, (3i is the switching rate from state 1 to 2, while (32 determines 
the transition rate 2^1. Note that the generalization for the renovation processes with arbitrary probability densities 
for switching times is straightforward. An important feature of the present analysis is an observation of the influence 
of the migration-proliferation dichotomy on the overall invasion rate of cancer cells. In what follows we show how the 
overall propagation rate u depends on the parameters (3k. 

We consider the growing tumor spheroid consisting of the tumor core with a high density of cells and the outer 
invasive zone where the cell density is much smaller. To describe the cancer cells of the two phenotypes we introduce 
the density of the cells of migrating phenotype, ni (t , x) , and the density of the cells of proliferating phenotype, n2 (i , x) . 
The balance equations for ni(t, x) and n2(t, x) are 

ni{t,x) = ni(0,x)\E'(t)e"'^i* -h / / ni{t - s,x- z)^{s,z)e~^^'dzds 

Jo jR-i 

+/32 / n2(t-s,x)*(s)e-^i*ds, (2) 
Jo 

n2(i,x) = n2(0,x)e-^=*-F / f {ni{t- s,^),n2{t- s,yi))e-^^'ds 

Jo 

+(3i [ m{t- s,yi)e-'^^'ds, (3) 
Jo 

where $(s,z) is the joint probability density function of making a jump z in the time interval s to s -|- ds, and 
denotes the integration is over d-dimensional space. The one dimensional case {d = 1) was considered in [7]. 

Cell migration (random jumps) involves a receptor-mediated adhesion to matrix proteins, matrix degradation by 
proteases, detachment from adhesion sites, active invasion into "new" intercellular space formed by degradation, etc. 
It would be extremely difficult to build up a rigorous deterministic model for this process. Since these factors are 
too many, we believe that a good alternative to such a model is a random walk with memory eff'ects. The active 
mechanism of migration of tumor cells involves small random jumps and delay time between jumps. The latter might 
be of the same order as the proliferation time. This dynamics is obviously random and its distribution is given by the 
probability density function (pdf) ^(s) : 

i^{s)= [ $(s,z)dz, (4) 

Jr'I 

where $(s, z) is the joint pdf. 

Equation (2) is the conservation law for cells of type 1 at time t at position x. The first term on the right hand side 
ni(0,x)\I/(t)e^'^i' represents cells of type 1 that stay up to time t at position x such that no jump occurred, and no 
switch took place. This term involves the function ^'(f) 



^{t) = 1 - / tlj{s)ds 
Jo 



(5) 
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which is the probabiUty that a cell of type 1 makes no jump until time t . Note that the exponential factor 

e-/3-=* = i_ f v{Tk)ds, k = l,2 
Jo 

is the probability that cells of phcnotypcs k do not switch until time t. The independence of the random jumps and 
switching gives us the probabihty \l/(t)e^'^i* while the first factor ni(0, x) is the initial density of cells of type 1 at x. 
The second term 

/ / ni(i - s,x - z)$(s,z)e~''^''(iz<is 

Jo JR-i 

gives us the number of cells of type 1 arriving at x up to time t. We assiimc the following random mechanism of 
migration: the cell of type 1 at time t — s sX position x — z waits a random time s before jumping a distance z at 
position x and remains a cell of type 1. The last term 

132 I n2{t- s,x)^{s)e-'^"'ds 
Jo 

represents the number of cells of type 2 that switch to the cell of type 1 up to time t and remain the cells of type 1 
(the factor e^'^i"). It also takes into account the fact that if transition 2—^1 happens at time t — s, then no jump 
takes place during the remaining time s (the factor 5'(s)). 

Equation (3) describes the balance of cells of proliferating phenotypc (no jumps). The first term on the right hand 
side, n2(0,x)e~'^2*, is the density of cells of type 2 that stay up to time t at position x such that no switch 2^1 
takes place. The second term on the right hand side 

t 

f {ni{t - s,x),n2{t - s,x)) e-f^^^ds 

is the proliferation rate for cell of type 2, which occurs providing that no switch takes place up to time t. The last 
term 

Pi [ nj_{t- s,x)e-^^^ds (6) 
Jo 

gives the number of cell of type 1 switching to the state 2 over the time interval (0, t). 

It is well known that the CTRW modelling is a standard technique for studying anomalous diffusion [8, 9]. We 
employ this technique to take into account subdiffusion that leads to slow motility of cancer cells in the invasive zone. 
In this paper each random step of a cancer cell is characterized by a waiting time s and a jump z which are distributed 
according to the joint pdf ^{s, z). This pdf can be written in a decoupled form 

$(s,z) = ^(s)p(z), (7) 

where ip{s) is waiting time pdf and p(z) is the pdf of cell jumps. This form corresponds to the case when the 
random waiting time and the individual displacement are independent. The subdiffusion regime occurs when the 
mean waiting time < t >= Tip{T)dT is infinite and the spherically symmetrical pdf /5(|x|) = p{r) has a finite 
variance = J r^p(r)dr < oo, where r is the radius of the spheroid. If the asymptotic behavior for the waiting-time 
density tp {t) for large t is t~^~^ with < C < 1) the mean waiting time < t > is infinite and the mean-square 
displacement a'^t'' corresponds to subdiffusion [8, 9]. When < t > in finite, there is normal diffusion: the mean-square 
displacement is Dt, where = cr^ /6 < t > for the three-dimensional case (d = 3). Superdiffusion takes place when 
the variance is infinite. Note that in many of the superdiffusion realization cases, the decoupling assumption of 
Eq. (7) can be inappropriate [17]. In what follows we consider only two regimes: normal diffusion and subdiffusion. 



B. Integro-differential equations 



The interesting feature of the balance equations (2) and (3) with $(s,z) = ^(s)p(z) is that they can be rewritten 
as a system of integro-differential equations: 

^= / a{t-s)[ [ni{s,x-z)-ni{s,x)]p{z)dzds- I3ini+f32n2, (8) 
Ci Jo JR'' 
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= f{ni,n2) + Pmi - I32n2, (9) 

where the memory kernel a(t) has to be determined. Let us derive these equations from (2) and (3) by using the 
Laplace transform for il>{t), and the Fourier transform for p(x) 



(10) 



Jo Jr<^ 
and the Fourier-Laplace (F-L) transform for the densities nk{t,x) 

fik{H,'k) = TC[nk{t,^)\= [ [ nfe(t,x)e-^*+'''-^didx, k=l,2. (11) 

Jr<' Jo 

Equation (2) with $(s,z) = tp{s)p{z) in the F-L space reads 

MH,-k) = ni(Q,k) ^ ~ + + h,{H, k)p(k)Vi(g + Pi) 

(12) 

To perform the F-L transform in Eq. (11) we use the standard convolution property 



rii(if,k)p(k)v;(if) = j 1^' 



ni{t — s,x — z)p{z)t{j{s)dzds 



dtdy 



Rearranging Eq. (12) and introducing the 'memory' kernel a{t) in term of its Laplace transform: 

H + Pi)i>{H + l3i) 
(l-Vi(if + /3i)) 



aiH) = i^i±MM±^, (13) 



we obtain 

Hhi{H,-k) - ni(0,k) = ^i(ff,k)a(ff)(p(k) - 1) + /32^2(i?,k) - /3i^i(F,k) . (14) 

Applying the F-L transform inversion to Eq. (14), we obtain the integro-differential equation (8). To find the F-L 
transform of Eq. (3), we denote the nonlinear proliferation term by Z{t,x) = f (ni(t, x), n2(i, x)). Its F-L transform 
is 

2{H,k)=£J^[Z{t,x)] . (15) 

We have from Eq. (3) 

MH,]^) = MO,^)jj^^ + kH,^)jj^^+PMH,^)jj^^ , (16) 

where Z{H,'k)/{H + P2) = JOJ^ /q Z{t - s,x)e-f^^''ds. Rearranging Eq. (16) in the following form 

Hfi2{H, k) - n2(0, k) = Z{H, k) + piMH, k) - P2n2{H, k) , 
and applying the F-L inversion and using Eq. (15), we obtain Eq. (9). 



C. Probability density function for cell jumps 

Now we are in a position to discuss different approximations for the probability density function for cell jumps p(z). 
Of course this function is not symmetrical in general. The cells of the migrating phenotype are biased to migrate away 

from the tumor spheroid core. The reasons for this asymmetrical creeping arc the non-uniform nutrient concentration 
(chemotaxis) , the gradient of cell adhesion sites (haptotaxis) , etc. Experimental observations suggest that cell jumps 
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are controlled by adhesion of tumor cells to extracellular matrix and jump lengths are very small [2]. Therefore p(z) is 
a rapidly decaying function for large |z|. In other words, the density of tumor cells varies on the scales that are much 
larger than the typical jump length. Thus one can use the Taylor series in Eq. (2) with $(s, z) = ip{s)p{2.) expanding 
ni(t — s,x— z) in z and truncate the series at the 2nd moment. This truncation for rapidly decaying function p{z) is 
a well defined procedure, since the higher moments become progressively smaller [18]. We have 

f , N / X , / X dn\ 1 d'^rii 

j ni{t - 5,x - z)p(z)dz = ni{t - s,x)- < Zi > — + ^ < ^i^o > dx dx- """ ' ^ ' 

where the Einstein rule for summation over repeated indices i and j is implied, and angular brackets denote averaging 
with respect to p{z) : 

< Zi >= I Zip{z)dz , < ZiZj >= / ZiZjp{z)dz . (18) 
Substitution of Eq. (17) into Eq. (2) with the decouple property $(s,z) = ?/)(s)p(z) yields 

ni(i,x) = ni(0,x)*(i)e-^i*+ j ni{t - s,x)i;{s)e-'^''ds- < Zi > j ^'tl;{s)e-'^''ds 

Jo Jo 

+1 < ZiZj> f ^^^{s)e-P''ds + l32 f n2{t-s,x)^{s)e-^''ds. (19) 



Note that the third term on the right hand side of this equation reflects a bias of random walk in the direction < z > . 
In fact, this equation involves the first two moments for random jumps: < Zi> and < ZiZj > . It can be rewritten as 
the integro-differential equation 

<Zi> a{t - s)-^ds = -< ZiZj > a{t - s) Q^^^^ ds - pim + /32n2 . (20) 

If the cell jumps are normally distributed then the characteristic function of p(z) is 

p(k) = exp ^iuiki - ^(Tijkikj^ , (21) 

where the summation convection is implied for the repeated index. The positive definite matrix aij can be written in 
terms of the first two moments 

aij =< ZiZj> — < Zi >< Zj> . (22) 

The probability density function p(z) is 



(23) 



where (cr ^)^^. is an inverse matrix. If wc assume that there is no bias, < Zj>=0, and < ZiZj >= for « ^ j, and 
< zf >= = ^. Then Eq. (19) takes the form 

ni(i,x) = ni(0,x)*(t)e-'^i* + / ni{t - s,x)il){s)e-^'^ds 

Jo 

+— Ani{t-s,x)il){s)e-'^''ds + p2 n2{t - s,x)^{s)e-'^''ds . (24) 
2« Jo Jo 



Prom the last equation one obtains integro-differential equation for m in d dimension 

v2 /•* 



dni a 
dt ^ub ./Q 



= 2dj - s)Ani{s,x)ds - pmi + /Jana . (25) 
Note that the one- dimensional case { d = 1) was analyzed in [7]. 
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D. Memory kernel and waiting time probability density function 

The formula 

gives us the relationship between the transport memory kernel a{t) in (20) and the waiting-time pdf 'tp{t) in terms of 
their Laplace transforms. It should be emphasized that it is impossible to find an explicit expression for the memory 
kernel a{t) for arbitrary choices of the waiting-time pdf tp{t). However, we are concerned with the rate of the spreading 
of tumor cells. In what follows we show that this rate depends on the Laplace transform a{H) rather than a{t). That 
is why the formula (26) is so important for our analysis. It follows from (26) that the transport kernel a{t) depends 
on the parameter /3i. This means that we can not separate the transport process and random switching in general. 
This phenomenon has been discussed recently in the literature on anomalous transport with reactions [19, 20]. 
Let us consider three typical distributions for the waiting-time pdf ip(t). 
(i) Exponential distribution. The random waiting time is exponentially distributed if it has a density 



ip{t) = Ae"^* . (27) 



The Laplace transform for this distribution is 



m) = I \e-^'e-"'dt = (28) 



and 



(l-^/;(^f + /3l)) 



therefore a{t) = X5{t). In this case the kernel a{t) is independent of Thus we have a classical system of convection- 
diffusion-reaction equations 

dni dn\ ^ 9^ni ^ _ 

'W'^'^'d^^ - Am + 132712 , (30) 



-^ = /(nin2)+/3ini-/32n2, (31) 

with the diffusion tensor Dij = A < ZiZj > /2 and the velocity v =A < z >. 

(ii) Gamma distribution. The waiting-time pdf V(i) corresponds to the family of gamma distributions with 
parameters m and A: 

\m4.ra — l^ — Xt 



Then il){H) = j and 



For example, if m = 2 



and the memory kernel is 



(g + /3i)A'" 

"^^^= (A + g + /?,)--A"^ - (''^ 



a{t) = A2e-(2A+/3i)t . (35) 
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The main result here is that the transport memory kernel depends on the parameter The integro-differential 
equation for cells of migratory phenotype takes the form 

^ + ..A /* e-^''+^^>^ds = A,A /* e-('^+^^>p^ds-P^m + 02n, . (36) 
ot J„ dxi Jo dxidxj 

The integro-differential Eq. (36) can be rewritten as the hyperbolic reaction-transport equation, and corresponding 
travelling wave solutions can be found as in [21, 22] (see also [14]). 

(Hi) Power law waiting time distribution,. The power law il^it) ^ 1/(1 + t/rY'^'' with < 7 < 1 is used in many 
applications [9] . Here we use r which is (in general case) not equal to 1/ A to stress the fractional property of cell 
dynamics. It is more convenient to use its Laplace transform 

= TIW- ^''^ 

Then 

^^^^ ^ + + y ^ iH±MZ. . (38) 

(l-^(if + /3i)) 

III. CANCER SPREADING RATE 

The overall rate u at which cancer cells spread is usually defined as the velocity of the experimentally detectable 
tumor front. In the generic Fisher equation setting the propagation rate is u = 2\/DU, where D is the diffusion 
coefficient and U is the proliferation rate [18]. The speed of this front is determined by the processes taking place 
at the leading edge of the cells' profile. In this paper we have a system of equations (2) and (3) and we define the 
overall spreading rate as the speed of the travelling wave solution of this system. For front-like initial conditions, the 
fronts for both densities ni and 712 quickly achieve the stationary forms that propagate with a constant rate u. The 
main purpose here is to find the dependence of this propagation rate on the statistical characteristics of the random 
switching process, (3i and /32. two moments for random jumps: {zi) and {ziZj) and waiting time distribution V'(i)- We 
use the logistic growth for cell proliferation 

/(m, na) = Un2 (1 - (ni + n2)/K) , (39) 

where U is the cell proliferation rate and K is the carrying capacity of the environment. We assume that the initial 
tumor spheroid of radius R has the following distribution of cells 

where positive constant Ai and A2 represent the stable equilibrium points of the densities ni and n2- They can be 
found from two equations Ai + A2 = K and PiAi = (32 A2 : 

I32K , /3iK , , 

We assume that the characteristic length scale for the tumor front is much smaller than the radius of the initial tumor 
spheroid. We also assume that the bias acts in the radial direction such that (z) = {r)er. These assumptions allow 
us to consider the propagation of the effective plane front in the radial direction, neglecting all curvature efi^ects. We 
expect that the long time development leads to the propagation of travelling fronts of permanent forms: ni (r — ut) 
and 712 (r — ut) , where the rate u is common to both densities ni and n2. 
The balance equations for densities ni and 712 are of the form 

ni(i,r) = ni(0,r)*(t)e-'^i*-K / ni{t - s,r)ip{s)e-'^''ds- < r > f ^i}{s)e-^^'ds 

Jo Jo 

+^ ^V'(s)e-^^^ds + 132 n2{t - s, r)^{s)e-^^'ds , (42) 
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n2{t,r) = n2{0,r)e-'^^* + U f n2{t - s,r){l - {ni{t - s,r) +n2{t - s,r)) / K)e-^^^ds 

Jo 

/ ni{t- s,r)e-'^^^ds. (43) 
Jo 

This system of equations is a starting point for the analysis of plane front propagation in a radial direction. 

A. Hyperbolic scaling and Hamilton-Jacobi equation 



The objective here is to find the rate u without resolving the shape of the travelling waves [12, 23]. For this purpose 
we use a hyperbolic scaling r r/e, t t/e and the rescaled density n| {t,r) = Uk {t/e,r/e) (see Appendix A). We 
write the density n| {t, r) in the exponential form 



nl {ti r) = Ak exp 



{t,r) 



A; = 1,2 



(44) 



where the non-negative function {t, r) describes the logarithmic asymptotic of both densities and plays a very 
important role. It follows from (44) that as long as the function 



G{t,r) = lim {t,r) 



£->0 



(45) 



is positive, the rescaled density n\{t,r) ^ as e ^ 0. We may argue that the equation G{t,r{t)) = gives 
us the spreading front position r {t) in the long-time and large-distance limit [12]. Substitution of the exponential 
transformation (44) into the equations for the rescaled densities n|(t, r) and taking the limit £ ^ yield two equations 
for Ai and A2. These equations have a non- trivial solution when the corresponding determinant is equal to zero (see 
Appendix A). It gives the following equation for G {t,r) : 



dG 



dG 



esr^^(s)e-/5i«ds 

J Jo 

poo POO 

■/3i/32 / e^"*(s)e-^i"rfs x / e^'e'^^'ds = 0. 
Jo Jo 



1 - U 



(46) 



In terms of the Laplace transform ip{H) = C[xl}{t)], Eq. (46) can be rewritten as a generalized Hamilton-Jacobi 
equation 



dG fdG 



^(-f + A) 



/3i/32(l-V;(-^ + /3i)) 



Note that inferring Eq. (47), we do not make any assumptions regarding waiting time pdf ip{s)- 



(47) 



B. Wavefront velocity 

Let us introduce the Hamiltonian function H and the generalized momentum p 

dG 



dr 



Then Hamilton-Jacobi equation (47) takes the form of the quadratic equation: 



<r > p + 



1 



1 



/3i/32(l-V'(g + /3i)) 
{H + /3i) {H + P2- U) 



2d V(-H" + /3i) 

This equation is very important because it allows us to find the spreading rate u 

dH H 



+ 1 = 0. 



u = 



dp p{H) • 



(48) 



(49) 



(50) 
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(51) 



We may equivalently write u = miriH | | , so u = , where H can be found from equation 

dp ^ p{H) 
dH H 

Let us illustrate this formula by using the classical Fisher equation 

dn ^d^n , 
^ = ^^ + Ml-n) 

for which the Hamiltonian is if = Dp'^/2 + U. Using this expression, we obtain 

,(„,= (H^_^2)"\ ,5.) 

Prom Eqs. (51) and (52) we obtain H = Dp^{H) = 2U, and therefore, the spreading rate for the Fisher equation is 
Up = H/p{H) = 2{DU)^ . This is the classical propagation speed. 

In what follows we consider a case when the mean jump length in the radial direction is zero, < r >= 0. If the 
random waiting time is exponentially distributed (27): tp{t) = Ae~^*, then the equation for the migratory cells is 

^ = D^-f3^m+f32n2. (53) 
The momentum p{H) can be found from (49) 

2 (g + /?l) JiJ2 

^ = — D DiH + p2-uy (''^ 

If we assume that /?i = P2, we can find from (50) that H = U and (54) p = (U/D)^^^, and H = U. There- 

1 /2 

fore, the spreading rate is uq = (UD) ' which is half of the classical Fishcr-KPP (Fisher-Kolmogorov-Petrovskii- 

Piskunov) propagation speed up. This result shows that the propagation rate is independent of the random migration- 

1/2 

proliferation switching for (ii — jd2- When /3i ^ j32 one can find the ratio of the propagation rate u and Uq = iUD) ' 
as 

u\ ' ^ 4(g + /32 - ^7)'^ \{H + (32 - U){H + p,) - p,p2] 
uoj [^H + p2-U)^+l3il32f 

where H is determined by Eq. (51). For the fixed values of /3i and U, the wavefront propagation rate versus /32//3i 
is depicted in Fig. 1. 

For the power law distribution ■)/'(t) ~ {r/tY^^ with < 7 < 1, the mean waiting time is divergent: < t >= 00. 
This assumption alone leads to the temporal fractional differential operator and corresponding anomalous diffusion 
equation [9]. The mean squared displacement for mobile cells is 

where = /2dT^ . 

Let us find the overall propagation of cancer cells as a result of interaction of subdiffusion (56), logistic proliferation 
(39), and random migration-proliferation switching (1). For the Laplace transform 'ii){H) = (1 -|- [Ht)^) ^ , the 
momentum p{H) can be found from (49): 

2 ^ (g + /3i)^ _ /3i/32(g + /?i)^-^ . . 

^ D^{H + l32-U)' ^ ' 

This formula together with (50) allows us to find the overall propagation rate of tumor cells u-y for the subdiffusion 
case. The case 7=1 corresponds to normal diffusion. One can find from (50), (54) and (57) the ratio of the anomalous 
propagation rate and the normal rate u determined by (55): 

^ = (H^T + (3iT)^ . (58) 
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Figure 1: (Color online) Propagation speed {u/uo)^ vs /32//3i- The values of Pi/U = [3, 4.5, 6.5] correspond to plots (1), (2), 
and (3), respectively. The insert corresponds to a solution of Eqs. (57) and (58) for the same values of (3i/U and 7 — 0.7. 



Since the "microscopic" time t is much smaller than the characteristic "cell proliferation" time and switching 
time and ~ J7, we conclude that Ht + Pit < 1. This condition of H ^ U is also confirmed by numerical 
solutions of Eqs. (51), (54), (57), and (58) (see insert in Fig. 1). It follows from (58) that the ratio u^/u increases 
up to 1 with 7 in the interval < 7 < 1. This means that normal diffusion leads to overestimation of the overall 
cancer spreading. Note that the advantage of balance Eqs. (2) and (3) is that they are related to a "mesoscopic" 
description of migratory cancer cells, and give us the statistical meaning of the phenomenological reaction-transport 
equation (20). 



IV. REACTION-TRANSPORT EQUATIONS 



The influence of the migration-proliferation dichotomy on the overall propagation rate is an important factor in 
glioma development. The Markovian switching between two phenotypes described by Eq. (1) can be generalized for 
the case when memory effects are taken into account. The system of integro-differential equations (8) and (9) takes 
the form 

a{t — s) / [ni{s,x — z) — ni{s,x.)] p{z)dzds 

I [Ai2(t - s)n2(s,x) - ^i(i - s)ni(s,x)]ds , (59) 





= /(ni,«2) - y [^2(^ - 5)712(5, x) - ^i(t " s)ni(s,x)]ds, (60) 

where /i^ (t) is the memory kernel for non-Markovian switching. Combining Eqs. (59) and (60) one finds that a total 
density n = ni + n2 obeys the equation 



^ = [ a{t~s) [ [ni(s,x - z) - ni(s,x)] p(z)dzds + /(ni,n2)- 
Jo Jr-^ 



(61) 



This equation does not restrict any possible random transitions between migration and proliferation phenotypes. 
Moreover, it can be a starting point of the glioma modelling in the framework of the differential equations. It can be 
rewritten in terms of the total density alone, if we introduce the probabilities pj such that ni = pin and n2 = P2n. By 
using the logistic growth for cell proliferation /(rii, 712) = /(^2) = Un2 (1 — n2/K) and rescaling p2n — > n, we obtain 



— = pi / a{t — s) / [n{s,x — z) — n{s,x)] p{z)dzds + Up2n {1 — n/ K) . 
Jo Jr-^ 



(62) 



11 



Let us find these probabilities for Markovian switehing (1). In fact there are four characteristic times in our model: 
proliferation time {U~^ for logistic growth), the transport time < t >= Jg tip{t)dt (averaging waiting time), and two 
switching times (3^^ and If we assume that both switching times are small compared to the growth time 
and transport time < t >, the "fast" switching process can be averaged. The "fast" local dynamics of densities ni 
and 712 governed by the equations 

dni a , o 
^ = -/3im+/32n2, ^ 



/3ini - p2n2 . 



(63) 



The solution for any x is 



ni(t) 



n2(i) = 



/32 



Pi +02 
01+02 



111 



(0) 



02 



n2 (0) 



01+02 

01 



01+02. 



-{01+02)t 



-(/3i+/32)t 



(64) 
(65) 



For the large intermediate time T such that 0^^ 0^^ <^T <^U ^, we have a local equilibrium, that is, ni = 



01+02 

If we consider now the transport and proliferation, it is clear that the total number of cancer cells n 



and n2 = 

splits locally to -^^q^n of migrating phenotype and -^^^n of proliferating phenotype. So 



01 



ni(i,x) = 



02 



01+0: 



-n{t,x), n2(t,x) 



01 



01+0: 



-n{t,x) . 



(66) 



This means that we have only one variable n{t, x) for which we can formulate a balance equation considering the 
transport for ni(t, x) and proliferation for n2{t,x). The probabilities are 



Pi 



02 



01+02' 



P2 



01 



01+02 



(67) 



In this limiting case, the model can be formulated in terms of the linear balance equation for the total number of 
cancer cells per unit volume n{t, x) 



n{t, x) 



02 



01+02 



Jo Jr'' 



01 /"* 

s,-x. — z)tp{s)p{z)dzds + — —U / n(t — s,x)ds. 

01 + 02 Jo 



(68) 



This reaction-transport equation can be also used to study the wavefront propagation in the framework of the 
Hamiltonian-Jacobi approach. 



V. CONCLUSION 



We developed a probabilistic approach for a migration-proliferation dichotomy in the spreading of tumor cells in 
the invasive zone. We derived the balance equations for densities of cancer cells of two phenotypes. In the migratory 
state the cells randomly move but there is no cell proliferation, while in the proliferating state the cancer cells do 
not migrate and only proliferation takes place. We took into account random switching between cell proliferation 
and migration by using a two-state Markov chain. The transport of tumor cells is formulated in terms of the CTRW 
with an arbitrary waiting time distribution, while proliferation is modeled by a non-linear function of both densities. 
We found the overall rate of tumor cell invasion for both normal diffusion and subdifFusion. The advantage of our 
probabilistic approach is that it allows us to take into account anomalous (subdiffusive) transport within the general 
scheme of migration, proliferation, and phenotype switching. We showed the equivalence of balance equations to 
a system of integro-differential equations involving memory effects for the transport of mobile cells. By using a 
hyperbolic scaling and Hamilton-Jacobi formalism we derived formulae for the overall spreading rate of cancer cells. 
We showed that the memory effects (subdiffusion) leads to a decrease in propagation rate compared to a standard 
diffusion approximation for transport. 
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Appendix A 

Rescaling of Eqs. (42) and (43), we obtain 

n\{t,r) = nf(0,r)*(t)e-''^*/^+ / nl{t - es,r)ip{s)e-'^''ds 

Jo 

+/32 / nl{t - es,r)^{s)e-^^^ds, 
Jo 



(69) 



ft/e 

n|(i,r) = n|(0, Oe"'^^*/^ + [/ / n|(i - £s,r)(l - (nf + n|) /ii')e-'^="ds 

Jo 



rt/s 

/ nf(i — £s,r)£ 
Jo 



-P^^ds. 



(70) 



Substitution of the exponential transformation n|(f, r) = A/j exp ^— ^-^^^ into these equations and accounting 
initial conditions yields 



Ax 



pt/e 

= Ai exp 
Jo 



tlj{s)e-'^''ds 



-e < r>Ai exp ( ) / — exp ^ — V(s)e '^^^ds 



dr 



+- 



2d 



■ exp 



+P2A2 I exp 
/o 



G-(/, ;■) - G"^(/ - j.s. /•) 



*(s)e"'^i"ds, 



(71) 



t/e 



+l3iAx I 
Jo 



A2 = UA2 I exp 

rt/e 



G^{t,r)-G%t-£s,r) 



exp 



G^{t,r)- G^{t- es,r) 



Ax + A2 ( G' 



-^■"'ds . 



-^^'ds 



Taking the limit e — > we have 



Ai = Ai e^'^{s)e-^''ds + Ai<r>^ / ew«^(s)e-/3i^rfs 
Jo c/r Jo 



(72) 



(73) 



A2 = UA2 / e^'e-^^'ds + /JiAi / e^'e-^^'ds 
Jo Jo 



Then Eqs. (73) and (74) can be rewritten as a system of linear equations for Ai and A2 

Ai 



( dG fdGV\ r ea, , . 

/•OO 

-A2/32 / e W^^f(s)e-/3i^rfs = , 
Jo 



(74) 



(75) 



/o 



1-U [' e^'e-^^'ds 
Jo 



0. 
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This system has a non-trivial solution when the corresponding determinant is equal to zero: 



dG 



1- l+<->^+^u7 



dG 



■L 



l-U e^'e-^^'ds 



pOO POO 

-/3i/32 / e^'^{s)e-l^^'dsx e^'e'^^'ds = Q . 
Jo Jo 



(77) 
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